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ABSTRACT 



Force-free magnetic fields are important in many astrophysical settings. De- 
termining the properties of such force-free fields - especially smoothness and 
stability properties - is crucial to understanding many key phenomena in astro- 
physical plasmas, for example energy release processes that heat the plasma and 
lead to dynamic or explosive events. Here we report on a serious limitation on 
the computation of force-free fields that has the potential to invalidate the results 
produced by numerical force-free field solvers even for cases in which they appear 
to converge (at fixed grid resolution) to an equilibrium magnetic field. In the 
present work we discuss this problem within the context of a Lagrangian relax- 
ation scheme that conserves magnetic flux and V • B identically. Error estimates 
are introduced to assess the quality of the calculated equilibrium. We go on to 
present an algorithm, based on re-writing the curl operation via Stokes' theorem, 
for calculating the current which holds great promise for improving dramatically 
the accuracy of the Lagrangian relaxation procedure. 

Subject headings: magnetic fields — methods: numerical — stars: coronae 
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Introduction 



Force-free magnetic fields, B, satisfy J x B = 0, or equivalently 

(V X B) X B = 0. ( 
Calculation of such force-free fields is of importance in many astrophysical settings, for 



example accretion disks around various objects (e.g. Frank et al. 2002 Uzdensky et al. 



2002 


), neutron stars ( 


McKinney 


2006), pulsars (Mestel 


1973 


), magnetic clouds ( 


Burlaga 



1988), and solar and stellar coronae (e.g. Anzer 1968). 



A particular application in solar physics is the controversial 'topological dissipation' 



model proposed by Parker (1972). The assertion of this model is that if an equilibrium 
magnetic field is perturbed by arbitrary motions at a line-tied boundary, then the 
subsequent field cannot relax to a smooth force-free equilibrium. Rather, the equilibrium 
must contain tangential discontinuities - corresponding to current sheets. Doubt has been 
cast upon the model however, as a number of authors have demonstrated the existence of 
smooth solutions in the scenario posed ( van Ballegooijen||1985 Zweibel fc Li||1987 Longcope 



& Strauss 1994; Craig & Sneyd 2005). The question as to whether current sheets form 
spontaneously in the coronal magnetic field is key to understanding the so-called coronal 
heating problem. This is just one example which demonstrates that determining both the 
structure and stability of force-free magnetic fields is of fundamental importance. 

There are different approaches that one may take when searching for force-free magnetic 
fields. One method, often used when modelling the solar corona, is to solve a boundary value 



problem (Amari et al. (1997), and see Schrijver et al. (2006) for a comparison of numerical 



schemes). The force- free field is reconstructed from boundary data, provided for example 
by a vector magnetogram. An alternative approach is to begin with an initial magnetic field 
that is not force-free and to perform a relaxation procedure. This is the natural approach 
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if one wants to investigate the properties of particular magnetic topologies. As long as the 
relaxation procedure can be guaranteed to be ideal, then the topology will be conserved 
during the relaxation. 

One powerful computational approach for investigating the properties of force-free 
fields is to employ an ideal Lagrangian relaxation scheme. Such schemes exploit the 
property that under ideal MHD the vector B/p evolves according to the equation 

D B _ /B 

where D/Dt is the material derivative, p the plasma density and v the plasma velocity. 
This is of exactly the same form as the evolution equation of a line element (5x in a flow 
(see, e.g. 



(2) 



Moffatt (1978)), and thus a Lagrangian description facilitates a relaxation that is. 



by construction, ideal. These schemes can be used to investigate the structure and (ideal 
MHD) stability of force-free fields. The latter is guaranteed by the iterative convergence 
of the scheme provided that the resolution is sufficient. The primary variables that the 
numerical scheme dynamically updates are the locations of the mesh points, with the 
quantities B and J being calculated via matrix products involving the initial magnetic field 
and derivatives of the mapping that describes the mesh deformation. An artificial frictional 



term is included in the equation of motion (see also Chodura & Schlueter (1981)) which 



guarantees a monotonic decrease of the energy. Two implementations of this method are 



described in Craig & Sneyd (1986) and Longbottom et al. (1998). The method has been 



used extensively to investigate the stability and equilibrium properties of various different 



magnetic configurations, such as the kink instability of magnetic fiux tubes (Craig & Sneyd 



1990), line-tied collapse of 2D and 3D magnetic null points (Craig & Litvinenko 2005 



Pontin & Craig 2005) and the Parker problem (Longbottom et al. 1998 Craig & Sneyd 



2005). 



In the following section we describe a test problem that illustrates one major difficulty 
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in the computation of force-free fields, in the context of the Lagrangian relaxation scheme 
outlined above. In Section [3] we present two possible extensions of the numerical scheme. 
In Section |4] we describe our results, and in Section [5] we present our conclusions. 



2. The problem 
2.1. Outline of the problem 



In a numerical relaxation experiment using braided initial fields (Wilmot-Smith et al 



2009) we came across an inconsistency of the resulting numerical force-free state, which 
is best explained with the help of the following example. Consider a magnetic field 
obtained from the homogenous field by a simple twisting deformation as shown in Fig. [l](a). 
Obviously an ideal relaxation towards a force-free state must end again in a homogenous 
state (J = 0). During this process the Lagrangian relaxation leads to a deformation of 
the initial computational mesh which exactly cancels the initial deformation applied to 
the homogenous field. This is a well defined setup in which we know exactly the initial 



and final states. We now employ the implicit (ADI) relaxation scheme detailed by Craig 



& Sneyd (1986) to relax our twisted field to a force-free equilibrium. The magnetic field 
is line-tied on all boundaries (B ■ n = on x and y boundaries). The J x B force as 
calculated by the numerical scheme decreases monotonically to an arbitrarily small value 
(e.g. 10~^-10~^), giving the appearance that the scheme converges (in an iterative sense) 
to a force- free equilibrium (to any desired accuracy, down to machine precision). However, 
when plotting a, the force-free proportionallity factor, along a field line it shows variations 
which are by orders of magnitude higher than would be expected from |J x B| < 10^^. It is 
this inconsistency that we investigate in what follows. As we will discuss the convergence of 
the numerical scheme in what follows, it is worth emphasising here the distinction between 
iterative convergence (at fixed resolution N) and real convergence, i.e. convergence towards 
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a 'correct' solution as the resolution N becomes sufficiently large. 



2.2. Analysis 



In order to investigate the source of the inconsistency described in the previous section, 
we consider the test problem outlined there. Specifically, we begin with an initially uniform 
magnetic field (B = boz), and superimpose two regions of toroidal field, centred on the 
2;- axis at ±zq, with exactly the same functional form, but of opposite signs: 



regions of twisted field, which are of opposite sign, should exactly cancel one another under 
an ideal relaxation, approaching the uniform field (with J = 0) as the equilibrium. Note 
that 1 01 is the maximum turning angle of field lines around the z-a.xis. 

One of the great advantages of an ideal Lagrangian relaxation is that it is possible to 
extract the paths of the magnetic field lines of the final state if one knows them in the 
initial state, simply by interpolating over the mesh displacement. Calculating the field lines 
in this way, no error is accumulated by integrating along B. Given knowledge of the field 
line paths, one can test the quality of the force-free approximation by plotting a along field 
lines. For a force- free field 



and a should be constant along field lines since taking the divergence of the above yields 




(3) 



with 60 = 1, ctr = V^, ciz = and 0i = vr, 02 = — tt, Li = —4, L2 = 4. We refer to this 
field in the following as T2. The field T2 (see Fig. [l|^a)) is constructed such that the two 



V X B = aB 



(4) 



B ■ Va = 0. 



(5) 
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(b) 

Fig. 1. — (a) Sample field lines for the field T2, given by Eq. ([s]). (b) Mesh in the z = 
plane for the test problem with artificially imposed deformation, with ip = n and resolution 
81^. 
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We begin by defining the variable a*, motivated by Eq. Q, as 

«* = ^- (6) 

We find that for the magnetic field calculated by the relaxation scheme, the value of a* 
changes dramatically along field lines. Of course the relaxation gives a magnetic field for 
which J X B is not identically zero. So for a given value of J x B, what is the maximum 
possible variation in a* along a field line? 

Consider 

V- J = V-(J||B) + V- Jx = 5, 
say, where S is representative of the error in calculating J. Using Eq. (j6| to replace J|| gives 

B-Va* = -V ■J± + 6. (7) 

or 

da* V ■ J_L 5 

IT - -IbT W\- 

where / is a parameter along a magnetic field line with units of length. Now suppose that 
|J X B|/|Bp < e within our domain. This implies that |J^| < e |B|, so that 

|V.JJ<^'^' 



d ' 

where d is the length scale of variations perpendicular to the magnetic field. Then from 
Eq. dSl) 



da* 



dl 



Returning to our relaxation results, we have for example e = 10~^, with |B| ^ 1, 
d ^ \/2. However, we find that \da* / dl\.fnax ~ 0.02. The discrepancy between this figure 
and the value of e must come from the final term in Eq. ([o]). This has been checked by 
interpolating the data onto a rectangular mesh and approximating V ■ J using standard 
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finite differences. We find V ■ J ~ 0(10"^), and it tfierefore appears tliat tlie residual 
currents parallel to B are not relaxed because V ■ J 7^ 0. As demonstrated below, this error 
does however decrease as the resolution is increased (see Tables [T]-[4]) . 

It turns out that the appearance of the errors is related to the way in which J 
is calculated within the scheme, via a combination of 1st and 2nd derivatives of the 
deformation matrix. These derivatives are calculated via finite differences in the numerical 
scheme, and it is here that these discretisation errors arise. This is demonstrated below. 



2.3. Accuracy test: artificially imposed deformation 

To ascertain the source of the errors, we take our initial state T2 and instead of 
performing the relaxation procedure, we artificially apply a deformation to the mesh which 
we can write down as a closed form expression, and moreover for which we can obtain the 
derivatives of the mesh displacement, and thus the resultant B and J fields, as closed form 
expressions. Motivated by the results of the relaxation, we impose a similar rotational 
distortion of the mesh which acts to 'untwist' the field, via the transformation 

{x,y,z) — > {xcos9 — ysinO , ?/ cos 6* + x sin 6* , z) (10) 

where 

« = V,exp(-^^-g. (11) 

ip constant. We now apply this transformation to an initially rectangular mesh on which 
B is given by T2, and compare the numerical and exact values for each entry in the mesh 
deformation Jacobian, and each component of B and J. Results are shown for three 
different values of the parameter ip in Table [TJ 

It is clear that there are large errors in the current calculated by the numerical scheme. 
While errors in B and in each individual term in the mesh distortion Jacobian are much 
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N 


= 71 


= 37r/4 


ip = 7r/2 


21 


29.5 


15.3 


6.46 


670 


159 


37.8 










41 


7.22 


3.93 


1.72 


144 


42.4 


10.6 










61 


3.30 


1.80 


0.793 


64.7 


19.9 


4.99 










81 


1.86 


1.02 


0.451 


37.0 


11.3 


2.83 



Table 1: Errors in B and J for deformations with ip = it, 37r/4, 7r/2 (in Eq. (11)) using 2nd- 
order finite differences. A^'^ is the mesh resolution. In each case the upper number shows the 
maximum relative percentage error in the domain over all components of B, i.e. 100 x \Bi — 



Bi\max/\Bf\max, where B" is the exact value. The lower number is 100 x \ Ji- Ji\max/\Ji\max- 
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smaller, it turns out that the combination in which they are multiplied, summed, and 
divided to calculate J incurs large errors. The calculation has been meticulously checked 
such that we are certain that the error appears not due to mathematical or coding error, but 
rather due to an accumulation of numerical truncation errors in the process of calculating J 

. Note that the errors increase as the mesh distortion 
{■ip) increases, and decrease with resolution (A^). 



(via Eq. (2.10) in Craig fc Sneyd|[l986 



3. More sophisticated numerical schemes 

3.1. Higher-order derivatives 

Clearly the accuracy of the force-free approximation is impaired by the accuracy of curl 
operation performed in the numerical scheme. One way to increase the accuracy of spatial 
derivatives could be to use higher-order finite difference expressions. Existing versions of the 
scheme use conventional second-order centred-differences involving two nearest-neighbour 
(n.n.) values. If we expect smooth solutions (without grid-scale features) then increasing to 
fourth-order finite difference expressions (using 4 n.n.) is expected to increase the accuracy. 

Recall that in the frictional Lagrangian method fluid displacements are determined 
from an equation of the form 

where A and C are prescribed tensor and vector functions and summation over repeated 
(Greek) indices is assumed. Note that partial differentiation with respect to the 
background Cartesian cordinates {Xi, X2, X^) is indicated using the comma notation (i.e. 
d'f/d(3d^ = U,). 

The important point for us is that the method involves two spatial derivatives of the 
Lagrangian variables Xa- This reflects the fact that the Lorentz force is computed using first 
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and second order derivatives of the Lagrangian mesh. Now "diagonal derivatives" such as 
Xijj are relatively easy to compute using finite differences: they involve the point itself and 
two/four nearest neighbours depending on whether the scheme is second or fourth order. 
It is these derivatives that are handled implicitly (via tri-diagnonal and penta-diagonal 
implementations of the ADI method) to provide, formally at least, the unconditional 
stability of the numerical scheme. Note that computation of the mixed derivatives can be 
more complicated: terms such as Xjjfc involve sixteen terms in the fourth order scheme, as 
opposed to just four when the method is second order. However, irrespective of the formal 
accuracy, perhaps the main drawback of the scheme is that, unlike V ■ B, the numerical 
evaluation of V ■ J is not guaranteed to vanish identically. Thus the accuracy of the relaxed 
solution can be compromised by the presence of rogue currents especially in weak field 
regions where the mesh is highly distorted. The examples presented below show explicitly 
that this can restrict the convergence of the solution with resolution A^. 



3.2. A routine based on Stokes' theorem 



Here we present an algorithm for calculating the curl of a vector field (say J = V x B) 
on a non-uniform mesh. This algorithm gives promising results, as discussed below, and is 
based on re-writing the curl operation, via Stokes' theorem, as a line integral: 



B dr 



J • fi (iS* 



J„ dS 



(12) 



'C JU{C) JU(C) 

where the surface U{C), with unit normal n, has the closed curve C as its boundary. 



The idea is similar to that of Hyman & Shashkov (1997) who have applied such so-called 



'mimetic' numerical methods to solving Maxwell's equations (Hyman & Shashkov 1999). 
Our algorithm differs in some ways from theirs. 



Eq. (12) can be discretised as follows. Suppose that we want to calculate J at the 



Fig. 2. — Notation used for calculation of J via the Stokes-based routine. 

mesh point Xij^k- There are three mesh surfaces that intersect at this point. The first is 
the ith mesh level in the first index direction. Consider the circuit in this surface shown 
in Fig. |2| and let the nearest neighbour points to j ^ be denoted x/, x//, x///, x/y 
Then, defining dxi = x// — x/, dx2 = x/// — x/j, etc. and Bi = (B(x/) + B(x//))/2, 
B2 = (B(x//) + B(x///))/2 etc., we can approximate the left hand side of Eq. ([I2|) by 



J = Bi • dxi + B2 ■ dx2 + B3 ■ dx3 + B4 ■ dx4. (13) 

Furthermore, the area of the enclosed quadrilateral is 

1111 
A = - |dxi X dx2| + - |dx2 X dxsl + - |dx3 x dx4| + - |dx4 x dxi| (14) 

We can now define the direction perpendicular to this mesh surface as 

(•j^^ 1 f dxi X dx2 dx2 X dx3 ^ dx3 x dx4 ^ dx4 x dxi 



4 V|dxi X dx2| |dx2 X dx3| |dx3 x dx4| |dx4 x dxi| ' 



and (comparing Eqs. (12 14)) the current in this direction as 

J(^) = I /A. (16) 

In the same way, we can define JA and Jn perpendicular to the other two mesh surfaces, 
with normal vectors n*^^) and n^^\ that pass through Xij^i^. 
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Now, the Jn^ are projections of the current we require (denoted J^) in such a way that 

n(^)-J, = J(^), p= 1,2,3. (17) 
Denoting by TV the matrix whose rows are the row vectors n^^^ n^^^ and n^^\ and by J„ 



the vector with components Jn \ Jn \ Jn^ we can re-write Eq. (17) as M^s = Jn- Finally, 
we obtain the current at point ^ ^ via 

J, =Ar-ij„. (18) 

Note that M is always invertible assuming that the n^^^ are linearly independent, i.e. as long 
as the grid cells have non-zero volume. It is straightforward to verify that this procedure 
reduces to the standard 2nd-order centred difference expression in each direction for a 
rectangular (undeformed) mesh. 

Since this method makes use of Stokes' theorem, which is "topological" in the sense 
that is does not depend on a deformation of the loop we are integrating over, we expect the 



method to be more robust and accurate for deformed grids (see also Hyman & Shashkov 



1997). The algorithm has not yet been implemented in the full numerical scheme, as it 
requires a complete re-writing of the implicit time-stepping routine, and a simple explicit 
implementation turns out to be prohibitively slow to run at reasonable resolution. However, 
the algorithm is tested and used as a diagnostic in what follows. 



4. Comparison of methods 

We now return to the test problem with two equal and opposite twists centred on 
z = ±zo described above. 
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4.1. Relaxing towards a rectangular mesh 

First, to demonstrate tliat in principle the original (2nd-order) relaxation scheme is 
sound, we perform the relaxation 'experiment' in the following way. We begin with a 
uniform field (B = boz) on a uniform mesh, and then deform this mesh in such a way that 
the resulting magnetic field is T2 (with (pi = ivr). This configuration is then relaxed, to a 
level where the Lorentz force calculated by the numerical scheme is reduced to J x B = 10^^. 
The equilibrium field corresponding to T2 is the uniform field B = boz, and in this case 
(due to the frozen-in condition) the straight field should correspond to the rectangular 



mesh. (The choice |0| = vr is motivated by the braiding example discussed in Wilmot-Smith 



et al. (2009) since this is the minimum level of twist which yields a magnetic field whose 



field lines are truly braided in that case.) 

To diagnose the success of the relaxation, referring back to Eq. ([9]), we calculate 

taking d = \/2 (radius of twist regions) and Aa* and A/ as the maximum change in a* over 
a given length, which occurs along the central field line - the z-axis - by symmetry. This 
expression puts a true value on the 'quality' e* of the force-free approximation: for a given 
variation in a*, e* provides a lower bound for the maximum value of |J x B|/|Bp within 
our domain, for a current free of errors (i.e. setting 5 = in Eq. ([o])). We obtain a value 
of e* = 9.5 X 10^^ for resolution = 61, demonstrating that discretisation errors in the 
scheme are very small when the relaxed state has an approximately rectangular mesh. 



4.2. Accuracy test: artificially imposed deformation 



We now investigate the promise of the two extensions to the scheme described in 
the previous section. In our test case (T2) the topology of the field is simple, and the 
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equilibrium field known, permitting the above approach (i.e. relaxation towards a uniform 
mesh). However, to investigate magnetic fields with non-trivial topology we must approach 
the problem in a different way. We therefore return to the case where we begin with T2 on 
a rectangular mesh, setting 0j = ±7r. 

Performing the artificially imposed analytical ('untwisting') deformation instead of 
relaxation, we see that derivatives calculated with the two new methods (4th-order finite 
differences and the Stokes-based method) both give smaller maximum errors than with 
the original 2nd-order scheme (see Table [2]). For a less deformed mesh {tp = 7r/2), the 
4th-order finite differences perform better than the Stokes-based routine. However, for a 
more distorted mesh {ip = tt), the Stokes routine gives significantly lower errors than either 
of the finite-difference methods. It is particularly interesting to note that the errors for the 
Stokes-based method seem to scale relatively weakly with the mesh deformation, suggesting 
it to be a good choice for highly deformed meshes. Both new schemes, for reasonable 
resolution and levels of deformation (A^ > 41, ip = it) give errors that are an order of 
magnitude lower than the original (2nd-order) method. This suggests these methods are 
worth pursuing, so we go on to perform relaxation simulations using them. 



4.3. Relaxation with initially rectangular mesh 

We now leave the artificially imposed deformation, and relax T2 using both the 2nd- 
and 4th-order schemes. The relaxation is allowed to run until |J x B| as calculated by the 
relevant numerical scheme is reduced to 10^^. We compare this value with e* (defined by 



Eq. (19)), and also the maximum value of the Lorentz force obtained by calculating J via 



the Stokes-based routine, denoted J<j. 

The results for two levels of initial twist (0j = ±|, (pi = ±7r) are displayed in Tables [s] 



(a) ' (b) ^ 

Fig. 3. — Isosurfaces of |J| at 2/3 of maximum, for (a) an intermediate stage in the relaxation 
(J X B = 0.05) and (b) the final state (J x B = 10"^). 4th-order scheme, N = 81, (pi = ±n. 
Inset: view in xy-plane (i.e. from z > 10). 

and|4j A number of points are immediately clear. First, in no simulation do we approach 
the apparent value of e = 10~^. However, the use of fourth- rather than second-order 
finite differences improves the quality of the relaxed field by an order of magnitude in e* 
when 101 = it/2. Increasing the deformation in the final state (by increasing |0| in the 
initial state T2) has a strong adverse effect on the relaxation process. This is found to be 
because spurious (unphysical) current concentrations arise where none should reasonably 
be expected. Examining the corresponding mesh, we find that these 'false' current regions 
appear where the grid is most distorted - see Fig. [3| and compare with Fig. [l|b). The 
'current shards' shown in Fig. |3](b) actually intensify as the relaxation proceeds. We find 
that this is possible since V ■ J (approximated by interpolating J onto a uniform mesh) is 
not close to zero. As a result there is no 'return current' associated with these localised 
current regions, which might be expected to generate a Lorentz force that would act against 
the further intensification of the current shards (if they have no physical basis). In previous 
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1^ = 71 




tjj = 7r/2 


N 


2 n.n. 


4 n.n. 


Stokes 




2 n.n. 


4 n.n. 


Stokes 


















21 


670 


222 


22.1 




37.8 


12.2 


18.7 


41 


144 


19.4 


5.85 




10.6 


1.01 


4.71 


61 


64.7 


4.14 


2.69 




4.99 


0.551 


2.13 


81 


37.0 


1.61 


1.57 




2.83 


0.898 


1.20 



Table 2: Errors in J for deformations with ijj = n and tt/2 (in Eq. (11)) using 2nd- and 
4th-order finite differences (2 n.n./4 n.n., respectively) and the Stokes-based routine. A^^ is 
the mesh resolution. In each case the value shown is the maximum relative percentage error 
in the domain over all components of J, i.e. 100 x | Jj — Ji\max/\Ji\max- 



N 


2 n.n. 




4 n.n. 


e* 


J. xB/|Bp 




e* 


JsXB/|B|2 














21 


0.11 


0.053 




0.048 


0.063 


41 


0.054 


0.046 




0.0067 


0.018 


61 


0.032 


0.035 




0.0042 


0.0050 


81 


0.021 


0.027 




0.0015 


0.0019 



Table 3: Values of the force-free quality parameter e* and the Lorentz force calculated via 
the Stokes-based routine, x B/|Bp, for simulation runs with 2nd- and 4th-order finite 
differences (2 n.n. / 4 n.n.) and resolution A^^, to two significant figures. Twist parameter 
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studies using such codes (e.g. Craig &: Litvinenko||2005 Pontin fc Craig||2005 ), intensification 
of |J| as |J X B| decreased in time was associated with current singularities, so at first sight 
it appears that these current shards could naively be interpreted as 'current sheets', which 
would of course be unphysical. Note, however, that the most important signature of current 
singularity in previous studies was a (power-law) proportionality of the peak current with 
mesh resolution (for given e). We have found that in fact the current shards become less 
intense as is increased, so there is a clear distinction between the two phenomena. 



Finally, consider the values of x B/|Bp we find for the relaxed fields (Tables [3]-|4]). 
They are clearly of the same order as e* (note that e* based on J from the numerical scheme 
and e* based on are of the same order), and thus it seems that an implementation 
involving the Stokes-based routine has the capacity to yield a magnetic field that is much 
closer to being force-free (with lower e*). This is illustrated in Fig. |4| We consider the 
observed value of a* based on Eq. ^ (dashed line), compared with a maximum allowable 
value for a*. Since a* is anti-symmetric about z = 0, the maximum value allowed, based 
on Eq. (g, is 

and we take |B| = 1 and d = \/2 as before, and L = 20, the length of the domain. Then we 



N 



21 
41 
61 
81 



2 n.n. 



0.28 
0.16 
0.11 
0.074 



JsXB/|B| 



0.18 
0.17 
0.13 
0.11 



4 n.n. 



J. xB/|B| 



0.17 
0.071 
0.026 
0.021 



Table 4: As Table [3} with twist parameter (pi = ±7r. 



0.21 
0.13 
0.062 
0.023 



- 20 - 




Fig. 4. — Evolution of different a*'s through the relaxation with second order finite differ- 
ences with = 81 and 0j = ±7r/2. The dashed line is the observed value of a*, while the 



solid line is the maximum allowable a* defined via Eq. (20) with e given by J x B/|B| from 
the 2nd-order numerical scheme. The dot-dashed line is the maximum a* with e given by 
3s X B/|Bp. Inset: close-up of behaviour at early time. 

obtain the maximum possible a* by taking e to be the maximum value during the relaxation 
of J X B (solid line) or x B (dot-dashed). We see that very early in the relaxation the 
actual value of a* becomes greater than the maximum allowed by J x B from the numerical 
scheme (2nd-order). Moreover, the discrepancy grow steadily. However, a* always remains 
less than the maximum allowed by the Stokes-based method, implying that this may be a 
more sound method to calculate the current and resulting Lorentz force. 
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5. Conclusions 

Force-free magnetic fields are important in many astrophysical applications. 
Determining the properties of such force-free fields - especially smoothness and stability 
properties - is key to understanding energy release processes that heat the plasma and lead 
to dynamic events such as flares in the solar corona. We have investigated the properties 
of different relaxation procedures for determining force-free fields based on a Lagrangian 
mesh approach. These techniques have previously been shown to have many powerful and 
advantageous properties. Previous understanding was that such schemes would iteratively 
converge (i.e. J x B decreasing monotonically to a given level) up to a certain degree of 
mesh deformation. Beyond this level of mesh deformation the scheme no longer converges 
(J X B oscillates or grows), and it is this phenomenon that was thought to limit the method. 
However, we have shown above that even when the numerical scheme iteratively converges, 
the accuracy of the force-free approximation can become seriously compromised for even 
'moderate' mesh deformations. This error is an accumulation of numerical discretisation 
errors resulting from the calculation of J via combinations of 1st and 2nd derivatives of the 
mesh deformation Jacobian - which are calculated using finite differences. The result is 
that neither J = V x B nor subsequently V • J = are well satisfied. 

It was demonstrated that a result of the breaking of the solenoidal condition for J can 
be the development of spurious (unphysical) current structures. However, we note that these 
rogue currents do diminish with resolution {N), so when using these schemes this property 
should always be checked where possible. We expect that, as a result, if it were possible to 
systematically increase N indefinitely the rogue currents would eventually vanish. In other 
words, the real problem is that the iterative convergence (i.e. monotonic decrease of |J x B| 
with t) is not compromised by the rogue currents but the real convergence to a correct 
solution is severely impaired. 
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A force-free field is defined by V x B = aB. One key result of this equation is that a 
must be constant along magnetic field lines. We therefore argued that a correct diagnostic 
to measure the quality of a force-free apprioximation is the constancy of the parameter 



a* = J||/|B| along field lines. An appropriate normalisation is given in Eq. (19). The 
results of our investigations suggest that for Lagrangian schemes the |J x B| measure does 
not provide a good indicator of true convergence — -a better measure in the constancy of a* 
along B. We note that other authors have proposed measures other than the maximum 



of J X B for testing a force-free approximation - for example Wheatland et al. (2000) 
introduced the "mean current-weighted angle between J and B". However, calculation 
of this measures still relies upon the value of J x B in the numerical scheme, and in the 
present scenario we have shown that the errors arise not because J and B are not parallel, 
but because V • J 7^ 0. 

Since errors in the (Lagrangian) numerical scheme investigated here arise as the mesh 
becomes increasingly distorted, a natural choice is to begin with a non-equilibrium field 
on a non-rectangular mesh, and relax towards a (perhaps approximately) rectangular one. 
However, this approach is not feasible if the field has complex topology. In the case of a 
braided field - which is of particular interest to the theory of the solar corona - we find that 



for our realisation of such a field (Wilmot-Smith et al. 2009) there is no escaping having at 



least a moderately distorted mesh in the final state. 

We proposed two possible extensions to the numerical method. The first was to 
increase the order of the finite differences used. It was found that for certain levels 
of deformation this can give an order of magnitude improvement in the quality of the 
force-free approximation obtained. It is therefore certainly a good approach to use in some 
circumstances. As the mesh became more and more highly deformed, the advantage of the 
scheme with 4th-order finite differences was lost for our test case T2. Furthermore, we 
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found that for relaxation of the braided field described in Wilmot-Smith et al. (2009) no 



appreciable improvement arose from using the 4th-order scheme. 

The other extension that we proposed to the scheme seems very promising. In Section 



3.2 we presented an algorithm for calculating the curl of a vector field on an arbitrary 
mesh, based on Stokes' theorem. For increasing levels of mesh deformation, this performed 
progressively better than the finite difference methods. What's more, in all of our tests 
the resultant Lorentz force x B had lower errors than that calculated by the traditional 
finite difference. In order for a relaxation experiment to remain accurate as it proceeds, the 



maximum allowed value of a* based on J x B (see Eq. (20)) must always remain greater 
than the maximum observed value of a* . We found that this is the case for the Stokes-based 
a* down to at least an order of magnitude lower in x B than for the finite difference 
methods (see Fig. |4]). 

All of the above leads us to believe that the Stokes-based algorithm is a highly 
promising one for improving the accuracy of Lagrangian relaxation schemes. At present it 
has not been implemented (i.e. the code does not act to minimise x B) because this 
requires a complete re-writing of the implicit (ADI) time-stepping, and a simple explicit 
implementation turns out to be prohibitively computationally expensive. However, our 
intended next step in this investigation is to implement this scheme, either by introducing 
the Stokes-based current calculation as a correction term in the existing scheme or by 
employing a more sophisticated explicit time-stepping to reduce the computational expense 
to acceptable levels. We note that while the algorithm at present only uses two nearest 
neighbour points in each direction, it could be extended to include further line integrals as 
corrections to the present formula for in much the same way as is done by increasing the 
order of finite difference derivatives. 



The authors are grateful to A. Nordlund and K. Galsgaard for helpful discussions. 
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